Minor Project: Oocyte Alternative Splicing Analysis

1 Introduction to this Analysis

1.1 Part 1: Upstream Analysis

The upstream analysis is performed on the CRG cluster using the Bash command line. This phase involves:

1. Data Acquisition:
   • Downloading *.fastq.gz files from 3 different databases.
2. Alignment:
   • Aligning the sequences to the latest mouse genome release (mm10) using Bowtie2 through the Vast-Tools align function.
3. Event Detection:
   • Importing inclusion tables generated by Vast-Tools into RStudio. The tables cover the following splicing events:
     • Exons (EX)
     • Introns (IN)
     • Alternative 5’ Splice Sites (Alt5)
     • Alternative 3’ Splice Sites (Alt3)
     • Microexons (MIC)

---

1.2 Part 2: Statistical Analysis

1.2.1 Objective

To identify splicing events that are differentially spliced in Oocytes treated with PladB with two differentes doses, or modified to be FMN2 DKO or SPIRE1/2 DKO.

1.2.2 Steps

1. Statistical Testing:
   • Use the betAS package in R to perform:
     • False Discovery Rate (FDR)
     • F-statistic
     • Probability of Differential Splicing (Pdiff)
2. Filtering Criteria:
   • Retain events with:
     • FDR ≤ 0.1
     • Pdiff ≥ 0.90
3. Batch Effect Mitigation:
   • A final list is generated by intersecting the splicing events identified in all five databases. This ensures the events are associated exclusively with differentiation. All exon related events (C1,C2,C3,S,including microexons) will be collapsed into a single list, EX.

---

1.3 Importing Inclusion Data

pladb_data <- getDataset(pathTables = paste0(getwd(),"/inclusion_tables/pladb_INCLUSION_LEVELS_FULL-mm10.tab"), tool = "vast-tools")
pladb_events <- getEvents(pladb_data, tool = "vast-tools") # Extract alternative splicing events
pladb_events <- alternativeEvents(pladb_events, minPsi = 1, maxPsi = 99) # Filter events based on PSI thresholds
pladb_exons <- filterEvents(pladb_events, types = c("C1", "C2", "C3", "S", "MIC"), N = 3)
pladb_alt <- filterEvents(pladb_events, types = c("Alt5", "Alt3"), N = 3)
pladb_introns <- filterEvents(pladb_events, types = c("IR"), N = 3)


fmndko_data <- getDataset(pathTables = paste0(getwd(),"/inclusion_tables/fmndko_INCLUSION_LEVELS_FULL-mm10.tab"), tool = "vast-tools")
fmndko_events <- getEvents(fmndko_data, tool = "vast-tools") # Extract alternative splicing events
fmndko_events <- alternativeEvents(fmndko_events, minPsi = 1, maxPsi = 99) # Filter events based on PSI thresholds
fmndko_exons <- filterEvents(fmndko_events, types = c("C1", "C2", "C3", "S", "MIC"), N = 3)
fmndko_introns <- filterEvents(fmndko_events, types = c("IR"), N = 5)
fmndko_alt <- filterEvents(fmndko_events, types = c("Alt5", "Alt3"), N = 3)



spire_data <- getDataset(pathTables = paste0(getwd(),"/inclusion_tables/spiredko_INCLUSION_LEVELS_FULL-mm10.tab"), tool = "vast-tools")
spire_events <- getEvents(spire_data, tool = "vast-tools") # Extract alternative splicing events
spire_events <- alternativeEvents(spire_events, minPsi = 1, maxPsi = 99) # Filter events based on PSI thresholds
spire_exons <- filterEvents(spire_events, types = c("C1", "C2", "C3", "S", "MIC"), N = 3)
spire_introns <- filterEvents(spire_events, types = c("IR"), N = 5)
spire_alt <- filterEvents(spire_events, types = c("Alt5", "Alt3"), N = 3)

1.4 Metadata

1.4.1 PladB

# Load metadata file containing sample information
metadata_pladb <- read.csv(paste0(getwd(),"/metadata/metadata_pladb.csv"), sep = "\t")
metadata_pladb<-metadata_pladb[seq(1, nrow(metadata_pladb), by = 2),]
DT::datatable(metadata_pladb, options = list(pageLength = nrow(metadata_pladb), scrollX = TRUE))

1.4.2 FMN2DKO

# Load metadata file containing sample information
metadata_fmn2dko<-data.frame(samples=c("SRR6026682_SRR6026683_SRR6026684_merged","SRR6026685_SRR6026686_SRR6026687_merged","SRR6026688_SRR6026689_SRR6026690_merged","SRR6026691_SRR6026692_SRR6026693_merged"),condition=c("control","control","fmndko","fmndko"))
DT::datatable(metadata_fmn2dko, options = list(pageLength = nrow(metadata_fmn2dko), scrollX = TRUE))

1.4.3 Spire

# Load metadata file containing sample information
metadata_spire<-data.frame(samples=c("Oocytes_FG_Spire12_Cont_a","Oocytes_FG_Spire12_Cont_b","Oocytes_FG_Spire12_Cont_c","Oocytes_FG_Spire12_DKO_a","Oocytes_FG_Spire12_DKO_b","Oocytes_FG_Spire12_DKO_c"), condition=rep(c("control", "spiredko"), each=3))
DT::datatable(metadata_spire, options = list(pageLength = nrow(metadata_spire), scrollX = TRUE))

---

2 Splicing & Sample Quality Checking

Splicing inclusion plots look like a U when looking at the exons, as normally most of them are spliced out completely (dPSI=0) or constitutively spliced in (dPSI=100), but some of them are in between and those are the most interesting for the analysis. On the other hand, introns are ussually not included (that is their definiiton), but in a few cases they are included. Thus, the intron plot looks like the exponential distribution.

2.1 PSI Distribution Plots

2.1.1 PladB

# Create the first plot for exon inclusion levels
bigPicturePlotExons <- bigPicturePlot(table = pladb_exons$PSI)

# Customize the first plot for professional reporting
plot1 <- bigPicturePlotExons +
  ggtitle("Exon Inclusion Levels") +
  xlab("Samples") +
  ylab("PSI Values") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    panel.background = element_rect(fill = "lightblue"),  # Set background color
    panel.grid.major = element_line(color = "gray90", size = 0.5),
    panel.grid.minor = element_line(color = "gray95", size = 0.25)
  )

# Create a second plot for a hypothetical dataframe
# Replace 'hypothetical_data$PSI' with your actual dataframe and column
bigPictureIntrons <- bigPicturePlot(table = pladb_introns$PSI)

# Customize the second plot
plot2 <- bigPictureIntrons +
  ggtitle("Introns Inclusion Levels") +
  xlab("Samples") +
  ylab("PSI Values") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    panel.background = element_rect(fill = "lightblue"),  # Set background color
    panel.grid.major = element_line(color = "gray90", size = 0.5),
    panel.grid.minor = element_line(color = "gray95", size = 0.25)
  )

# Combine both plots into one figure using patchwork
plot1 / plot2

2.1.2 Spire

# Create the first plot for exon inclusion levels
bigPicturePlotExons <- bigPicturePlot(table = spire_exons$PSI)

# Customize the first plot for professional reporting
plot1 <- bigPicturePlotExons +
  ggtitle("Exon Inclusion Levels") +
  xlab("Samples") +
  ylab("PSI Values") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    panel.background = element_rect(fill = "lightblue"),  # Set background color
    panel.grid.major = element_line(color = "gray90", size = 0.5),
    panel.grid.minor = element_line(color = "gray95", size = 0.25)
  )

# Create a second plot for a hypothetical dataframe
# Replace 'hypothetical_data$PSI' with your actual dataframe and column
bigPictureIntrons <- bigPicturePlot(table = spire_introns$PSI)

# Customize the second plot
plot2 <- bigPictureIntrons +
  ggtitle("Introns Inclusion Levels") +
  xlab("Samples") +
  ylab("PSI Values") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    panel.background = element_rect(fill = "lightblue"),  # Set background color
    panel.grid.major = element_line(color = "gray90", size = 0.5),
    panel.grid.minor = element_line(color = "gray95", size = 0.25)
  )

# Combine both plots into one figure using patchwork
plot1 / plot2

2.1.3 FMN2DKO

# Create the first plot for exon inclusion levels
bigPicturePlotExons <- bigPicturePlot(table = fmndko_exons$PSI)

# Customize the first plot for professional reporting
plot1 <- bigPicturePlotExons +
  ggtitle("Exon Inclusion Levels") +
  xlab("Samples") +
  ylab("PSI Values") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    panel.background = element_rect(fill = "lightblue"),  # Set background color
    panel.grid.major = element_line(color = "gray90", size = 0.5),
    panel.grid.minor = element_line(color = "gray95", size = 0.25)
  )

# Create a second plot for a hypothetical dataframe
# Replace 'hypothetical_data$PSI' with your actual dataframe and column
bigPictureIntrons <- bigPicturePlot(table = fmndko_introns$PSI)

# Customize the second plot
plot2 <- bigPictureIntrons +
  ggtitle("Introns Inclusion Levels") +
  xlab("Samples") +
  ylab("PSI Values") +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
    axis.text = element_text(size = 10),
    axis.title = element_text(size = 12),
    panel.background = element_rect(fill = "lightblue"),  # Set background color
    panel.grid.major = element_line(color = "gray90", size = 0.5),
    panel.grid.minor = element_line(color = "gray95", size = 0.25)
  )

# Combine both plots into one figure using patchwork
plot1 / plot2

2.2 Total Splicing Events

As it can be observed below, the dominic dataset has the highest number of splicing events, what can be explained by the higher number of reads (90M in myotubes and 140M in myoblasts).

# Adjusted version of the script
splicing_events <- tibble(
  event_type = names(spire_events$EventsPerType),
  pladb = pladb_events$EventsPerType,
  spire = spire_events$EventsPerType,
  fmndko = fmndko_events$EventsPerType
)

splicing_events_long <- splicing_events %>%
  pivot_longer(cols = c(pladb, spire, fmndko),
               names_to = "study", values_to = "count") %>%
  mutate(event_type = ifelse(event_type %in% c("C1", "C2", "C3", "ANN", "S", "MIC"), "EX", event_type)) %>%
  group_by(study, event_type) %>%
  summarize(count = sum(count, na.rm = TRUE), .groups = "drop")

# Proportion bar plot
plot_proportion <- ggplot(splicing_events_long, aes(fill = event_type, y = count, x = study)) +
  geom_bar(position = "fill", stat = "identity") +
  labs(
    title = "<b style='font-size:18px;'>Proportion of Splicing Events by Study</b>",
    x = NULL,
    y = "Proportion of Events",
    fill = "Event Type"
  ) +
  theme_minimal(base_family = font) +
  theme(
    axis.text.x = element_blank(),
    axis.ticks.x = element_blank(),
    plot.title = element_markdown(),
    panel.background = element_rect(fill = "gray98"),
    panel.grid = element_line(color = "gray90", size = 0.5)
  ) +
  scale_fill_paletteer_d("MetBrewer::Hokusai3")

splicing_events_long_sum<-splicing_events_long %>%
  dplyr::group_by(study) %>%
  summarize(count = sum(count))


# Absolute counts plot
plot_absolute <- ggplot(splicing_events_long_sum, aes(y = count, x = study)) +
  geom_point(color = "#C70039", size = 3) +
  geom_text(aes(label = count), vjust = -1, size = 4, check_overlap = TRUE) +
  labs(
    x = "Study",
    y = "Number of Events",
    caption = paste0("Created by AG on ", Sys.Date()), size=3
  ) +
  theme_minimal(base_family = font) +
  theme(
    legend.position = "none",
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.text.x = element_text( size=15),
    panel.background = element_rect(fill = "white"),
    panel.grid.major = element_line(color = "gray80", size = 0.5),
    panel.grid.minor = element_line(color = "gray95", size = 0.3)
  ) +
  scale_y_continuous(expand = expansion(mult = c(0.1, 0.2)))


# Combine plots
combined_plot <- plot_proportion / plot_absolute +
  plot_layout(heights = c(1, 0.5)) &
  theme(plot.margin = margin(5,5, 0, 5))

# Display combined plot
combined_plot

2.3 PCA Plot

2.3.1 PladB

# Subset and scale data
pca_pladb <- pladb_data[, c("EVENT",pladb_events$Samples)] %>%
  na.omit()
pca_spire <- spire_data[,  c("EVENT",spire_events$Samples)]
pca_fmndko <- fmndko_data[,  c("EVENT",fmndko_events$Samples)]
rownames(pca_pladb)<- pca_pladb$EVENT
pca_pladb<-t(pca_pladb[,-1])


pca_result_pladb <- prcomp(pca_pladb)  # Perform PCA

# Prepare PCA results for plotting
pca_data <- as.data.frame(pca_result_pladb$x)
pca_data$Sample <- rownames(pca_data)
pca_data$condition<-metadata_pladb$Description

# Create PCA plot
pca_plot <- ggplot(pca_data, aes(x = PC1, y = PC2, color = as.factor(condition), )) +
  geom_point(size = 6) +
  labs(
    title = "<b style='font-size:18px;'>PCA of Samples (PSI Data)</b>",
    x = paste("PC1 (", round(100 * summary(pca_result_pladb)$importance[2, 1], 1), "%)", sep = ""),
    y = paste("PC2 (", round(100 * summary(pca_result_pladb)$importance[2, 2], 1), "%)", sep = ""),
    caption = paste0("Created by AG on ", Sys.Date()), size=3
  ) +
  theme_minimal(base_family = "roboto") +
  theme(
    plot.title = element_markdown(),
    plot.title.position = "plot",
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    legend.position = "right",
    legend.title = element_blank(),
    panel.background = element_rect(fill = "gray98"),
    panel.grid = element_line(color = "gray90", size = 0.5),
    panel.grid.major = element_line(color = "gray90", size = 0.5),
    panel.grid.minor = element_blank()
  ) +
  scale_colour_paletteer_d("MetBrewer::Hokusai3") +
  coord_fixed()

pca_plot

2.3.2 SpireDKO

# Subset and scale data
pca_spire <- spire_data[, c("EVENT",spire_events$Samples)] %>%
  na.omit()

rownames(pca_spire)<- pca_spire$EVENT
pca_spire<-t(pca_spire[,-1])


pca_result_spire <- prcomp(pca_spire)  # Perform PCA

# Prepare PCA results for plotting
pca_data <- as.data.frame(pca_result_spire$x)
pca_data$Sample <- rownames(pca_data)
pca_data$condition<-metadata_spire$condition

# Create PCA plot
pca_plot <- ggplot(pca_data, aes(x = PC1, y = PC2, color = as.factor(condition), )) +
  geom_point(size = 6) +
  labs(
    title = "<b style='font-size:18px;'>PCA of Samples (PSI Data)</b>",
    x = paste("PC1 (", round(100 * summary(pca_result_pladb)$importance[2, 1], 1), "%)", sep = ""),
    y = paste("PC2 (", round(100 * summary(pca_result_pladb)$importance[2, 2], 1), "%)", sep = ""),
    caption = paste0("Created by AG on ", Sys.Date()), size=3
  ) +
  theme_minimal(base_family = "roboto") +
  theme(
    plot.title = element_markdown(),
    plot.title.position = "plot",
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    legend.position = "right",
    legend.title = element_blank(),
    panel.background = element_rect(fill = "gray98"),
    panel.grid = element_line(color = "gray90", size = 0.5),
    panel.grid.major = element_line(color = "gray90", size = 0.5),
    panel.grid.minor = element_blank()
  ) +
  scale_colour_paletteer_d("MetBrewer::Hokusai3") +
  coord_fixed()

pca_plot

2.3.3 FMN2DKO

# Subset and scale data
pca_fmndko <- fmndko_data[, c("EVENT",fmndko_events$Samples)] %>%
  na.omit()

rownames(pca_fmndko)<- pca_fmndko$EVENT
pca_fmndko<-t(pca_fmndko[,-1])


pca_result_fmndko <- prcomp(pca_fmndko)  # Perform PCA

# Prepare PCA results for plotting
pca_data <- as.data.frame(pca_result_fmndko$x)
pca_data$Sample <- rownames(pca_data)
pca_data$condition<-metadata_fmn2dko$condition

# Create PCA plot
pca_plot <- ggplot(pca_data, aes(x = PC1, y = PC2, color = as.factor(condition), )) +
  geom_point(size = 6) +
  labs(
    title = "<b style='font-size:18px;'>PCA of Samples (PSI Data)</b>",
    x = paste("PC1 (", round(100 * summary(pca_result_pladb)$importance[2, 1], 1), "%)", sep = ""),
    y = paste("PC2 (", round(100 * summary(pca_result_pladb)$importance[2, 2], 1), "%)", sep = ""),
    caption = paste0("Created by AG on ", Sys.Date()), size=3
  ) +
  theme_minimal(base_family = "roboto") +
  theme(
    plot.title = element_markdown(),
    plot.title.position = "plot",
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    legend.position = "right",
    legend.title = element_blank(),
    panel.background = element_rect(fill = "gray98"),
    panel.grid = element_line(color = "gray90", size = 0.5),
    panel.grid.major = element_line(color = "gray90", size = 0.5),
    panel.grid.minor = element_blank()
  ) +
  scale_colour_paletteer_d("MetBrewer::Hokusai3") +
  coord_fixed()

pca_plot

3 Differential Splicing Calculation

# Extract unique groups and sample IDs
# pladb

groupingVariable <- "Description"


groups_pladb <- unique(metadata_pladb[, groupingVariable])
samples_pladb <- metadata_pladb$fastq_files
samples_pladb <- paste0(
  "X",
  sub("\\.fastq\\.gz$", "", samples_pladb),
  "_merged"
)

# fmn2dko

metadata_fmn2dko<-data.frame(samples=c("SRR6026682_SRR6026683_SRR6026684_merged","SRR6026685_SRR6026686_SRR6026687_merged","SRR6026688_SRR6026689_SRR6026690_merged","SRR6026691_SRR6026692_SRR6026693_merged"),condition=c("control","control","fmndko","fmndko"))
groups_fmndko <- unique(metadata_fmn2dko[, "condition"])
samples_fmndko <- metadata_fmn2dko$samples

# spiredko

groups_spire <- unique(metadata_spire[, "condition"])
samples_spire <- metadata_spire$samples


random_colors=c("#D8D97AFF", "#95C36EFF", "#74C8C3FF", "#5A97C1FF", "#295384FF", "#0A2E57FF")

# Create group list with metadata
groupList_pladb <- lapply(1:length(groups_pladb), function(i) {
  list(
    name = groups_pladb[i],
    samples = samples_pladb[metadata_pladb[, groupingVariable] == groups_pladb[i]],
    color = random_colors[i]
  )
})
names(groupList_pladb) <- groups_pladb

# Create group list with metadata
groupList_fmndko <- lapply(1:length(groups_fmndko), function(i) {
  list(
    name = groups_fmndko[i],
    samples = samples_fmndko[metadata_fmn2dko[, "condition"] == groups_fmndko[i]],
    color = random_colors[i]
  )
})
names(groupList_fmndko) <- groups_fmndko

# Create group list with metadata
groupList_spire <- lapply(1:length(groups_spire), function(i) {
  list(
    name = groups_spire[i],
    samples = samples_spire[metadata_spire[, "condition"] == groups_spire[i]],
    color = random_colors[i]
  )
})
names(groupList_spire) <- groups_spire

# Define groups
groupA_pladb    <- "ControlFmn2+-"
groupC_pladb    <- "Fmn2+-_+1mMPlatB"
groupB_pladb <- "Fmn2+-_+10mMPlatB"

groupA_fmndko<-"control"
groupB_fmndko<-"fmndko"

groupA_spire<-"control"
groupB_spire<-"spiredko"

# Define samples inside each group
samplesA_pladb    <- groupList_pladb[[groupA_pladb]]$samples
samplesB_pladb    <- groupList_pladb[[groupB_pladb]]$samples
samplesC_pladb    <- groupList_pladb[[groupC_pladb]]$samples
colsGroupA_pladb    <- convertCols(pladb_exons$PSI, samplesA_pladb)
colsGroupB_pladb   <- convertCols(pladb_exons$PSI, samplesB_pladb)
colsGroupC_pladb   <- convertCols(pladb_exons$PSI, samplesC_pladb)

samplesA_fmndko <- groupList_fmndko[[groupA_fmndko]]$samples
samplesB_fmndko <- groupList_fmndko[[groupB_fmndko]]$samples
colsGroupA_fmndko <- convertCols(fmndko_exons$PSI, samplesA_fmndko)
colsGroupB_fmndko <- convertCols(fmndko_exons$PSI, samplesB_fmndko)

samplesA_spire <- groupList_spire[[groupA_spire]]$samples
samplesB_spire <- groupList_spire[[groupB_spire]]$samples
colsGroupA_spire <- convertCols(spire_exons$PSI, samplesA_spire)
colsGroupB_spire <- convertCols(spire_exons$PSI, samplesB_spire)


set.seed(42) #Setting seed for downstream simulations of the beta distribution

3.1 Exon Lists

3.1.1 Calculations

pladb_pdiff_exons <- prepareTableVolcano(
  psitable = pladb_exons$PSI,
  qualtable = pladb_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupB_pladb,
  labA = groupA_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_control_pdiff_exons <- prepareTableVolcano(
  psitable = pladb_exons$PSI,
  qualtable = pladb_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupC_pladb,
  labA = groupA_pladb,
  labB = groupC_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_10mm_pdiff_exons <- prepareTableVolcano(
  psitable = pladb_exons$PSI,
  qualtable = pladb_exons$Qual,
  npoints = 500,
  colsA = colsGroupC_pladb,
  colsB = colsGroupB_pladb,
  labA = groupC_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

# Prepare table for fmn2dko
fmndko_pdiff_exons <- prepareTableVolcano(
  psitable = fmndko_exons$PSI,
  qualtable = fmndko_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_fmndko,
  colsB = colsGroupB_fmndko,
  labA = groupA_fmndko,
  labB = groupB_fmndko,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

# Prepare table for spire
spire_pdiff_exons <- prepareTableVolcano(
  psitable = spire_exons$PSI,
  qualtable = spire_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_spire,
  colsB = colsGroupB_spire,
  labA = groupA_spire,
  labB = groupB_spire,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

pladb_fstat_exons <- prepareTableVolcanoFstat(
  psitable = pladb_exons$PSI,
  qualtable = pladb_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupB_pladb,
  labA = groupA_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_control_fstat_exons <- prepareTableVolcanoFstat(
  psitable = pladb_exons$PSI,
  qualtable = pladb_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupC_pladb,
  labA = groupA_pladb,
  labB = groupC_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_10mm_fstat_exons <- prepareTableVolcanoFstat(
  psitable = pladb_exons$PSI,
  qualtable = pladb_exons$Qual,
  npoints = 500,
  colsA = colsGroupC_pladb,
  colsB = colsGroupB_pladb,
  labA = groupC_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

# Prepare table for fmndko
fmndko_fstat_exons <- prepareTableVolcanoFstat(
  psitable = fmndko_exons$PSI,
  qualtable = fmndko_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_fmndko,
  colsB = colsGroupB_fmndko,
  labA = groupA_fmndko,
  labB = groupB_fmndko,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)


# Prepare table for spire
spire_fstat_exons <- prepareTableVolcanoFstat(
  psitable = spire_exons$PSI,
  qualtable = spire_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_spire,
  colsB = colsGroupB_spire,
  labA = groupA_spire,
  labB = groupB_spire,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

pladb_fdr_exons <- prepareTableVolcanoFDR(
  psitable = pladb_exons$PSI,
  qualtable = pladb_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupB_pladb,
  labA = groupA_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_control_fdr_exons <- prepareTableVolcanoFDR(
  psitable = pladb_exons$PSI,
  qualtable = pladb_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupC_pladb,
  labA = groupA_pladb,
  labB = groupC_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed=TRUE,
  CoverageWeight = FALSE)


pladb_1mm_10mm_fdr_exons <- prepareTableVolcanoFDR(
  psitable = pladb_exons$PSI,
  qualtable = pladb_exons$Qual,
  npoints = 500,
  colsA = colsGroupC_pladb,
  colsB = colsGroupB_pladb,
  labA = groupC_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed=TRUE,
  CoverageWeight = FALSE)


# Prepare table for fmndko
fmndko_fdr_exons <- prepareTableVolcanoFDR(
  psitable = fmndko_exons$PSI,
  qualtable = fmndko_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_fmndko,
  colsB = colsGroupB_fmndko,
  labA = groupA_fmndko,
  labB = groupB_fmndko,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed = TRUE,
  CoverageWeight = FALSE
)

# Prepare table for spire
spire_fdr_exons <- prepareTableVolcanoFDR(
  psitable = spire_exons$PSI,
  qualtable = spire_exons$Qual,
  npoints = 500,
  colsA = colsGroupA_spire,
  colsB = colsGroupB_spire,
  labA = groupA_spire,
  labB = groupB_spire,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed = TRUE,
  CoverageWeight = FALSE
)

3.1.2 FDR Volcano Plots

Volcano Plots of the FDR corrected exons for each study. In pink are the exons with a deltaPSI>=0.1.

3.1.2.1 PladB 10mM vs Control

plotVolcanoFDR(betasTable =filter(pladb_fdr_exons,!is.na(EVENT)),
                            labA = groupA_pladb,
                            labB = groupB_pladb,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +

theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.1.2.2 PladB 1mM vs Control

plotVolcanoFDR(betasTable =filter(pladb_1mm_control_fdr_exons,!is.na(EVENT)),
                            labA = groupA_pladb,
                            labB = groupC_pladb,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +

theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.1.2.3 PladB 10mM vs 1mM

plotVolcanoFDR(betasTable =filter(pladb_1mm_10mm_fdr_exons,!is.na(EVENT)),
                            labA = groupC_pladb,
                            labB = groupB_pladb,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +

theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.1.2.4 FMN2DKO

plotVolcanoFDR(betasTable =filter(fmndko_fdr_exons,!is.na(EVENT)),
                            labA = groupA_fmndko,
                            labB = groupB_fmndko,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +
theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.1.2.5 Spire2

plotVolcanoFDR(betasTable =filter(spire_fdr_exons,!is.na(EVENT)),
                            labA = groupA_spire,
                            labB = groupB_spire,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +
theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.1.3 Individual Combination Tables

These tables show the exons that show a FDR<=0.05 and PDiff (1-pvalue)>=0.95 in the pairwise comparison of each condition.

3.1.3.1 PladB 10mM vs control

pladb_pdiff_exons_tab <- pladb_pdiff_exons[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_pdiff_exons_tab) <- 1:nrow(pladb_pdiff_exons_tab)

pladb_fstat_exons_tab <- pladb_fstat_exons[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_fstat_exons_tab) <- 1:nrow(pladb_fstat_exons_tab)

pladb_fdr_exons_tab <- pladb_fdr_exons[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(pladb_fdr_exons_tab) <- 1:nrow(pladb_fdr_exons_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(pladb_pdiff_exons_tab$EVENT, pladb_fstat_exons_tab$EVENT, pladb_fdr_exons_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_exons_diff_tab <- pladb_pdiff_exons_tab %>% filter(EVENT %in% common_events)
filtered_exons_fstat_tab <- pladb_fstat_exons_tab %>% filter(EVENT %in% common_events)
filtered_exons_fdr_tab <- pladb_fdr_exons_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_exons_diff_tab <- filtered_exons_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_exons_fstat_tab <- filtered_exons_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_exons_fdr_tab <- filtered_exons_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
pladb_exon_combination_tab <- filtered_exons_diff_tab %>%
  inner_join(filtered_exons_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_exons_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Render DataTable with enhancements
datatable(
  pladb_exon_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(pladb_exon_combination_tab,
          file = paste0("pladb_exon_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.1.3.2 PladB 1mM vs Control

pladb_pdiff_exons_tab <- pladb_1mm_control_pdiff_exons[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_pdiff_exons_tab) <- 1:nrow(pladb_pdiff_exons_tab)

pladb_fstat_exons_tab <- pladb_1mm_control_fstat_exons[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_fstat_exons_tab) <- 1:nrow(pladb_fstat_exons_tab)

pladb_fdr_exons_tab <- pladb_1mm_control_fdr_exons[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(pladb_fdr_exons_tab) <- 1:nrow(pladb_fdr_exons_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(pladb_pdiff_exons_tab$EVENT, pladb_fstat_exons_tab$EVENT, pladb_fdr_exons_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_exons_diff_tab <- pladb_pdiff_exons_tab %>% filter(EVENT %in% common_events)
filtered_exons_fstat_tab <- pladb_fstat_exons_tab %>% filter(EVENT %in% common_events)
filtered_exons_fdr_tab <- pladb_fdr_exons_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_exons_diff_tab <- filtered_exons_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_exons_fstat_tab <- filtered_exons_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_exons_fdr_tab <- filtered_exons_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
pladb_1mm_control_exon_combination_tab <- filtered_exons_diff_tab %>%
  inner_join(filtered_exons_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_exons_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Find if events appear in the control to 10mM pairwise comparison
pladb_1mm_control_exon_combination_tab$appears_in_10mm_control <- pladb_1mm_control_exon_combination_tab$EVENT %in% pladb_exon_combination_tab$EVENT

# Render DataTable with enhancements
datatable(
  pladb_1mm_control_exon_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(pladb_1mm_control_exon_combination_tab,
          file = paste0("pladb_1mm_control_exon_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.1.3.3 PladB 10mM vs 1mM

pladb_pdiff_exons_tab <- pladb_1mm_10mm_pdiff_exons[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_pdiff_exons_tab) <- 1:nrow(pladb_pdiff_exons_tab)

pladb_fstat_exons_tab <- pladb_1mm_10mm_fstat_exons[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_fstat_exons_tab) <- 1:nrow(pladb_fstat_exons_tab)

pladb_fdr_exons_tab <- pladb_1mm_10mm_fdr_exons[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(pladb_fdr_exons_tab) <- 1:nrow(pladb_fdr_exons_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(pladb_pdiff_exons_tab$EVENT, pladb_fstat_exons_tab$EVENT, pladb_fdr_exons_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_exons_diff_tab <- pladb_pdiff_exons_tab %>% filter(EVENT %in% common_events)
filtered_exons_fstat_tab <- pladb_fstat_exons_tab %>% filter(EVENT %in% common_events)
filtered_exons_fdr_tab <- pladb_fdr_exons_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_exons_diff_tab <- filtered_exons_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_exons_fstat_tab <- filtered_exons_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_exons_fdr_tab <- filtered_exons_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
pladb_10mm_1mm_exon_combination_tab <- filtered_exons_diff_tab %>%
  inner_join(filtered_exons_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_exons_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Find if events appear in the control to 10mM pairwise comparison
pladb_10mm_1mm_exon_combination_tab$appears_in_10mm_control <- pladb_10mm_1mm_exon_combination_tab$EVENT %in% pladb_exon_combination_tab$EVENT

# Render DataTable with enhancements
datatable(
  pladb_10mm_1mm_exon_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(pladb_10mm_1mm_exon_combination_tab,
          file = paste0("pladb_10mm_1mm_exon_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.1.3.4 FMN2DKO

fmndko_pdiff_exons_tab <- fmndko_pdiff_exons[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(fmndko_pdiff_exons_tab) <- 1:nrow(fmndko_pdiff_exons_tab)

fmndko_fstat_exons_tab <- fmndko_fstat_exons[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(fmndko_fstat_exons_tab) <- 1:nrow(fmndko_fstat_exons_tab)

fmndko_fdr_exons_tab <- fmndko_fdr_exons[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(fmndko_fdr_exons_tab) <- 1:nrow(fmndko_fdr_exons_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(fmndko_pdiff_exons_tab$EVENT, fmndko_fstat_exons_tab$EVENT, fmndko_fdr_exons_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_exon_diff_tab <- fmndko_pdiff_exons_tab %>% filter(EVENT %in% common_events)
filtered_exon_fstat_tab <- fmndko_fstat_exons_tab %>% filter(EVENT %in% common_events)
filtered_exon_fdr_tab <- fmndko_fdr_exons_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_exon_diff_tab <- filtered_exon_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_exon_fstat_tab <- filtered_exon_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_exon_fdr_tab <- filtered_exon_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
fmndko_exon_combination_tab <- filtered_exon_diff_tab %>%
  inner_join(filtered_exon_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_exon_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Render DataTable with enhancements
datatable(
  fmndko_exon_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(fmndko_exon_combination_tab,
          file = paste0("fmndko_exon_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.1.3.5 Spire2

spire_pdiff_exons_tab <- spire_pdiff_exons[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(spire_pdiff_exons_tab) <- 1:nrow(spire_pdiff_exons_tab)

spire_fstat_exons_tab <- spire_fstat_exons[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(spire_fstat_exons_tab) <- 1:nrow(spire_fstat_exons_tab)

spire_fdr_exons_tab <- spire_fdr_exons[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(spire_fdr_exons_tab) <- 1:nrow(spire_fdr_exons_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(spire_pdiff_exons_tab$EVENT, spire_fstat_exons_tab$EVENT, spire_fdr_exons_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_exon_diff_tab <- spire_pdiff_exons_tab %>% filter(EVENT %in% common_events)
filtered_exon_fstat_tab <- spire_fstat_exons_tab %>% filter(EVENT %in% common_events)
filtered_exon_fdr_tab <- spire_fdr_exons_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_exon_diff_tab <- filtered_exon_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_exon_fstat_tab <- filtered_exon_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_exon_fdr_tab <- filtered_exon_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
spire_exon_combination_tab <- filtered_exon_diff_tab %>%
  inner_join(filtered_exon_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_exon_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Render DataTable with enhancements
datatable(
  spire_exon_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(spire_exon_combination_tab,
          file = paste0("spire_exon_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.2 Intron Lists

3.2.1 Calculations

pladb_pdiff_introns <- prepareTableVolcano(
  psitable = pladb_introns$PSI,
  qualtable = pladb_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupB_pladb,
  labA = groupA_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_control_pdiff_introns <- prepareTableVolcano(
  psitable = pladb_introns$PSI,
  qualtable = pladb_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupC_pladb,
  labA = groupA_pladb,
  labB = groupC_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_10mm_pdiff_introns <- prepareTableVolcano(
  psitable = pladb_introns$PSI,
  qualtable = pladb_introns$Qual,
  npoints = 500,
  colsA = colsGroupC_pladb,
  colsB = colsGroupB_pladb,
  labA = groupC_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

# Prepare table for Tao
fmndko_pdiff_introns <- prepareTableVolcano(
  psitable = fmndko_introns$PSI,
  qualtable = fmndko_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_fmndko,
  colsB = colsGroupB_fmndko,
  labA = groupA_fmndko,
  labB = groupB_fmndko,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

# Prepare table for spire
spire_pdiff_introns <- prepareTableVolcano(
  psitable = spire_introns$PSI,
  qualtable = spire_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_spire,
  colsB = colsGroupB_spire,
  labA = groupA_spire,
  labB = groupB_spire,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

pladb_fstat_introns <- prepareTableVolcanoFstat(
  psitable = pladb_introns$PSI,
  qualtable = pladb_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupB_pladb,
  labA = groupA_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_control_fstat_introns <- prepareTableVolcanoFstat(
  psitable = pladb_introns$PSI,
  qualtable = pladb_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupC_pladb,
  labA = groupA_pladb,
  labB = groupC_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_10mm_fstat_introns <- prepareTableVolcanoFstat(
  psitable = pladb_introns$PSI,
  qualtable = pladb_introns$Qual,
  npoints = 500,
  colsA = colsGroupC_pladb,
  colsB = colsGroupB_pladb,
  labA = groupC_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

# Prepare table for fmndko
fmndko_fstat_introns <- prepareTableVolcanoFstat(
  psitable = fmndko_introns$PSI,
  qualtable = fmndko_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_fmndko,
  colsB = colsGroupB_fmndko,
  labA = groupA_fmndko,
  labB = groupB_fmndko,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

# Prepare table for spire
spire_fstat_introns <- prepareTableVolcanoFstat(
  psitable = spire_introns$PSI,
  qualtable = spire_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_spire,
  colsB = colsGroupB_spire,
  labA = groupA_spire,
  labB = groupB_spire,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

pladb_fdr_introns <- prepareTableVolcanoFDR(
  psitable = pladb_introns$PSI,
  qualtable = pladb_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupB_pladb,
  labA = groupA_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_control_fdr_introns <- prepareTableVolcanoFDR(
  psitable = pladb_introns$PSI,
  qualtable = pladb_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupC_pladb,
  labA = groupA_pladb,
  labB = groupC_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed=TRUE,
  CoverageWeight = FALSE)


pladb_1mm_10mm_fdr_introns <- prepareTableVolcanoFDR(
  psitable = pladb_introns$PSI,
  qualtable = pladb_introns$Qual,
  npoints = 500,
  colsA = colsGroupC_pladb,
  colsB = colsGroupB_pladb,
  labA = groupC_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed=TRUE,
  CoverageWeight = FALSE)


# Prepare table for fmndko
fmndko_fdr_introns <- prepareTableVolcanoFDR(
  psitable = fmndko_introns$PSI,
  qualtable = fmndko_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_fmndko,
  colsB = colsGroupB_fmndko,
  labA = groupA_fmndko,
  labB = groupB_fmndko,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed = TRUE,
  CoverageWeight = FALSE
)

# Prepare table for spire
spire_fdr_introns <- prepareTableVolcanoFDR(
  psitable = spire_introns$PSI,
  qualtable = spire_introns$Qual,
  npoints = 500,
  colsA = colsGroupA_spire,
  colsB = colsGroupB_spire,
  labA = groupA_spire,
  labB = groupB_spire,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed = TRUE,
  CoverageWeight = FALSE
)

3.2.2 FDR Volcano Plots

Volcano Plots of the FDR corrected introns for each study. In pink are the exons with a deltaPSI>=0.1.

3.2.2.1 PladB 10mM vs Control

plotVolcanoFDR(betasTable =filter(pladb_fdr_introns,!is.na(EVENT)),
                            labA = groupA_pladb,
                            labB = groupB_pladb,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +

theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.2.2.2 PladB 1mM vs Control

plotVolcanoFDR(betasTable =filter(pladb_1mm_control_fdr_introns,!is.na(EVENT)),
                            labA = groupA_pladb,
                            labB = groupC_pladb,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +

theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.2.2.3 PladB 10mM vs 1mM

plotVolcanoFDR(betasTable =filter(pladb_1mm_10mm_fdr_introns,!is.na(EVENT)),
                            labA = groupC_pladb,
                            labB = groupB_pladb,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +

theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.2.2.4 FMN2DKO

plotVolcanoFDR(betasTable =filter(fmndko_fdr_introns,!is.na(EVENT)),
                            labA = groupA_fmndko,
                            labB = groupB_fmndko,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +
theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.2.2.5 Spire2

plotVolcanoFDR(betasTable =filter(spire_fdr_introns,!is.na(EVENT)),
                            labA = groupA_spire,
                            labB = groupB_spire,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +
theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.2.3 Individual Combination Tables

These tables show the introns that show a FDR<=0.05 and PDiff (1-pvalue)>=0.95 in the pairwise comparison of each condition.

3.2.3.1 PladB 10mM vs control

pladb_pdiff_introns_tab <- pladb_pdiff_introns[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_pdiff_introns_tab) <- 1:nrow(pladb_pdiff_introns_tab)

pladb_fstat_introns_tab <- pladb_fstat_introns[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_fstat_introns_tab) <- 1:nrow(pladb_fstat_introns_tab)

pladb_fdr_introns_tab <- pladb_fdr_introns[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(pladb_fdr_introns_tab) <- 1:nrow(pladb_fdr_introns_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(pladb_pdiff_introns_tab$EVENT, pladb_fstat_introns_tab$EVENT, pladb_fdr_introns_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_introns_diff_tab <- pladb_pdiff_introns_tab %>% filter(EVENT %in% common_events)
filtered_introns_fstat_tab <- pladb_fstat_introns_tab %>% filter(EVENT %in% common_events)
filtered_introns_fdr_tab <- pladb_fdr_introns_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_introns_diff_tab <- filtered_introns_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_introns_fstat_tab <- filtered_introns_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_introns_fdr_tab <- filtered_introns_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
pladb_intron_combination_tab <- filtered_introns_diff_tab %>%
  inner_join(filtered_introns_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_introns_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Render DataTable with enhancements
datatable(
  pladb_intron_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(pladb_intron_combination_tab,
          file = paste0("pladb_intron_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.2.3.2 PladB 1mM vs Control

pladb_pdiff_introns_tab <- pladb_1mm_control_pdiff_introns[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_pdiff_introns_tab) <- 1:nrow(pladb_pdiff_introns_tab)

pladb_fstat_introns_tab <- pladb_1mm_control_fstat_introns[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_fstat_introns_tab) <- 1:nrow(pladb_fstat_introns_tab)

pladb_fdr_introns_tab <- pladb_1mm_control_fdr_introns[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(pladb_fdr_introns_tab) <- 1:nrow(pladb_fdr_introns_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(pladb_pdiff_introns_tab$EVENT, pladb_fstat_introns_tab$EVENT, pladb_fdr_introns_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_introns_diff_tab <- pladb_pdiff_introns_tab %>% filter(EVENT %in% common_events)
filtered_introns_fstat_tab <- pladb_fstat_introns_tab %>% filter(EVENT %in% common_events)
filtered_introns_fdr_tab <- pladb_fdr_introns_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_introns_diff_tab <- filtered_introns_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_introns_fstat_tab <- filtered_introns_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_introns_fdr_tab <- filtered_introns_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
pladb_1mm_control_intron_combination_tab <- filtered_introns_diff_tab %>%
  inner_join(filtered_introns_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_introns_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Find if events appear in the control to 10mM pairwise comparison
pladb_1mm_control_intron_combination_tab$appears_in_10mm_control <- pladb_1mm_control_intron_combination_tab$EVENT %in% pladb_intron_combination_tab$EVENT

# Render DataTable with enhancements
datatable(
  pladb_1mm_control_intron_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(pladb_1mm_control_intron_combination_tab,
          file = paste0("pladb_1mm_control_intron_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.2.3.3 PladB 10mM vs 1mM

pladb_pdiff_introns_tab <- pladb_1mm_10mm_pdiff_introns[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_pdiff_introns_tab) <- 1:nrow(pladb_pdiff_introns_tab)

pladb_fstat_introns_tab <- pladb_1mm_10mm_fstat_introns[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_fstat_introns_tab) <- 1:nrow(pladb_fstat_introns_tab)

pladb_fdr_introns_tab <- pladb_1mm_10mm_fdr_introns[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(pladb_fdr_introns_tab) <- 1:nrow(pladb_fdr_introns_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(pladb_pdiff_introns_tab$EVENT, pladb_fstat_introns_tab$EVENT, pladb_fdr_introns_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_introns_diff_tab <- pladb_pdiff_introns_tab %>% filter(EVENT %in% common_events)
filtered_introns_fstat_tab <- pladb_fstat_introns_tab %>% filter(EVENT %in% common_events)
filtered_introns_fdr_tab <- pladb_fdr_introns_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_introns_diff_tab <- filtered_introns_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_introns_fstat_tab <- filtered_introns_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_introns_fdr_tab <- filtered_introns_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
pladb_10mm_1mm_intron_combination_tab <- filtered_introns_diff_tab %>%
  inner_join(filtered_introns_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_introns_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Find if events appear in the control to 10mM pairwise comparison
pladb_10mm_1mm_intron_combination_tab$appears_in_10mm_control <- pladb_10mm_1mm_intron_combination_tab$EVENT %in% pladb_intron_combination_tab$EVENT

# Render DataTable with enhancements
datatable(
  pladb_10mm_1mm_intron_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(pladb_10mm_1mm_intron_combination_tab,
          file = paste0("pladb_10mm_1mm_intron_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.2.3.4 FMN2DKO

fmndko_pdiff_introns_tab <- fmndko_pdiff_introns[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(fmndko_pdiff_introns_tab) <- 1:nrow(fmndko_pdiff_introns_tab)

fmndko_fstat_introns_tab <- fmndko_fstat_introns[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(fmndko_fstat_introns_tab) <- 1:nrow(fmndko_fstat_introns_tab)

fmndko_fdr_introns_tab <- fmndko_fdr_introns[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(fmndko_fdr_introns_tab) <- 1:nrow(fmndko_fdr_introns_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(fmndko_pdiff_introns_tab$EVENT, fmndko_fstat_introns_tab$EVENT, fmndko_fdr_introns_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_intron_diff_tab <- fmndko_pdiff_introns_tab %>% filter(EVENT %in% common_events)
filtered_intron_fstat_tab <- fmndko_fstat_introns_tab %>% filter(EVENT %in% common_events)
filtered_intron_fdr_tab <- fmndko_fdr_introns_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_intron_diff_tab <- filtered_intron_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_intron_fstat_tab <- filtered_intron_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_intron_fdr_tab <- filtered_intron_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
fmndko_intron_combination_tab <- filtered_intron_diff_tab %>%
  inner_join(filtered_intron_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_intron_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Render DataTable with enhancements
datatable(
  fmndko_intron_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(fmndko_intron_combination_tab,
          file = paste0("fmndko_intron_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.2.3.5 Spire2

spire_pdiff_introns_tab <- spire_pdiff_introns[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(spire_pdiff_introns_tab) <- 1:nrow(spire_pdiff_introns_tab)

spire_fstat_introns_tab <- spire_fstat_introns[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(spire_fstat_introns_tab) <- 1:nrow(spire_fstat_introns_tab)

spire_fdr_introns_tab <- spire_fdr_introns[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(spire_fdr_introns_tab) <- 1:nrow(spire_fdr_introns_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(spire_pdiff_introns_tab$EVENT, spire_fstat_introns_tab$EVENT, spire_fdr_introns_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_intron_diff_tab <- spire_pdiff_introns_tab %>% filter(EVENT %in% common_events)
filtered_intron_fstat_tab <- spire_fstat_introns_tab %>% filter(EVENT %in% common_events)
filtered_intron_fdr_tab <- spire_fdr_introns_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_intron_diff_tab <- filtered_intron_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_intron_fstat_tab <- filtered_intron_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_intron_fdr_tab <- filtered_intron_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
spire_intron_combination_tab <- filtered_intron_diff_tab %>%
  inner_join(filtered_intron_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_intron_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Render DataTable with enhancements
datatable(
  spire_intron_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(spire_intron_combination_tab,
          file = paste0("spire_intron_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.3 Alternative Splicing List

3.3.1 Calculations

pladb_pdiff_alt <- prepareTableVolcano(
  psitable = pladb_alt$PSI,
  qualtable = pladb_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupB_pladb,
  labA = groupA_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_control_pdiff_alt <- prepareTableVolcano(
  psitable = pladb_alt$PSI,
  qualtable = pladb_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupC_pladb,
  labA = groupA_pladb,
  labB = groupC_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_10mm_pdiff_alt <- prepareTableVolcano(
  psitable = pladb_alt$PSI,
  qualtable = pladb_alt$Qual,
  npoints = 500,
  colsA = colsGroupC_pladb,
  colsB = colsGroupB_pladb,
  labA = groupC_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

# Prepare table for Tao
fmndko_pdiff_alt <- prepareTableVolcano(
  psitable = fmndko_alt$PSI,
  qualtable = fmndko_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_fmndko,
  colsB = colsGroupB_fmndko,
  labA = groupA_fmndko,
  labB = groupB_fmndko,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

# Prepare table for spire
spire_pdiff_alt <- prepareTableVolcano(
  psitable = spire_alt$PSI,
  qualtable = spire_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_spire,
  colsB = colsGroupB_spire,
  labA = groupA_spire,
  labB = groupB_spire,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

pladb_fstat_alt <- prepareTableVolcanoFstat(
  psitable = pladb_alt$PSI,
  qualtable = pladb_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupB_pladb,
  labA = groupA_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_control_fstat_alt <- prepareTableVolcanoFstat(
  psitable = pladb_alt$PSI,
  qualtable = pladb_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupC_pladb,
  labA = groupA_pladb,
  labB = groupC_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_10mm_fstat_alt <- prepareTableVolcanoFstat(
  psitable = pladb_alt$PSI,
  qualtable = pladb_alt$Qual,
  npoints = 500,
  colsA = colsGroupC_pladb,
  colsB = colsGroupB_pladb,
  labA = groupC_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed=TRUE,
  CoverageWeight = FALSE)

# Prepare table for fmndko
fmndko_fstat_alt <- prepareTableVolcanoFstat(
  psitable = fmndko_alt$PSI,
  qualtable = fmndko_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_fmndko,
  colsB = colsGroupB_fmndko,
  labA = groupA_fmndko,
  labB = groupB_fmndko,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)


# Prepare table for spire
spire_fstat_alt <- prepareTableVolcanoFstat(
  psitable = spire_alt$PSI,
  qualtable = spire_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_spire,
  colsB = colsGroupB_spire,
  labA = groupA_spire,
  labB = groupB_spire,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  seed = TRUE,
  CoverageWeight = FALSE
)

pladb_fdr_alt <- prepareTableVolcanoFDR(
  psitable = pladb_alt$PSI,
  qualtable = pladb_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupB_pladb,
  labA = groupA_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed=TRUE,
  CoverageWeight = FALSE)

pladb_1mm_control_fdr_alt <- prepareTableVolcanoFDR(
  psitable = pladb_alt$PSI,
  qualtable = pladb_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_pladb,
  colsB = colsGroupC_pladb,
  labA = groupA_pladb,
  labB = groupC_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed=TRUE,
  CoverageWeight = FALSE)


pladb_1mm_10mm_fdr_alt <- prepareTableVolcanoFDR(
  psitable = pladb_alt$PSI,
  qualtable = pladb_alt$Qual,
  npoints = 500,
  colsA = colsGroupC_pladb,
  colsB = colsGroupB_pladb,
  labA = groupC_pladb,
  labB = groupB_pladb,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed=TRUE,
  CoverageWeight = FALSE)


# Prepare table for fmndko
fmndko_fdr_alt <- prepareTableVolcanoFDR(
  psitable = fmndko_alt$PSI,
  qualtable = fmndko_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_fmndko,
  colsB = colsGroupB_fmndko,
  labA = groupA_fmndko,
  labB = groupB_fmndko,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed = TRUE,
  CoverageWeight = FALSE
)

# Prepare table for spire
spire_fdr_alt <- prepareTableVolcanoFDR(
  psitable = spire_alt$PSI,
  qualtable = spire_alt$Qual,
  npoints = 500,
  colsA = colsGroupA_spire,
  colsB = colsGroupB_spire,
  labA = groupA_spire,
  labB = groupB_spire,
  basalColor = "#89C0AE",
  interestColor = "#E69A9C",
  maxDevTable = maxDevSimulationN100,
  nsim = 100,
  seed = TRUE,
  CoverageWeight = FALSE
)

3.3.2 FDR Volcano Plots

Volcano Plots of the FDR corrected Alt37Alt5 for each study. In pink are the exons with a deltaPSI>=0.1.

3.3.2.1 PladB 10mM vs Control

plotVolcanoFDR(betasTable =filter(pladb_fdr_alt,!is.na(EVENT)),
                            labA = groupA_pladb,
                            labB = groupB_pladb,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +

theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.3.2.2 PladB 1mM vs Control

plotVolcanoFDR(betasTable =filter(pladb_1mm_control_fdr_alt,!is.na(EVENT)),
                            labA = groupA_pladb,
                            labB = groupC_pladb,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +

theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.3.2.3 PladB 10mM vs 1mM

plotVolcanoFDR(betasTable =filter(pladb_1mm_10mm_fdr_alt,!is.na(EVENT)),
                            labA = groupC_pladb,
                            labB = groupB_pladb,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +

theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.3.2.4 FMN2DKO

plotVolcanoFDR(betasTable =filter(fmndko_fdr_alt,!is.na(EVENT)),
                            labA = groupA_fmndko,
                            labB = groupB_fmndko,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +
theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.3.2.5 Spire2

plotVolcanoFDR(betasTable =filter(spire_fdr_alt,!is.na(EVENT)),
                            labA = groupA_spire,
                            labB = groupB_spire,
                            basalColor = "#89C0AE",
                            interestColor = "#E69A9C") +
theme(
  plot.title = element_text(size = 10),  # Title font size
  axis.title = element_text(size = 10),  # Axis title font size
  axis.text = element_text(size = 10),   # Axis text font size
  legend.text = element_text(size = 10), # Legend text font size
  legend.title = element_text(size = 10) # Legend title font size
)

3.3.3 Individual Combination Tables

These tables show the Alt5 or Alt3 events that show a FDR<=0.05 and PDiff (1-pvalue)>=0.95 in the pairwise comparison of each condition.

3.3.3.1 PladB 10mM vs control

pladb_pdiff_alt_tab <- pladb_pdiff_alt[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_pdiff_alt_tab) <- 1:nrow(pladb_pdiff_alt_tab)

pladb_fstat_alt_tab <- pladb_fstat_alt[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_fstat_alt_tab) <- 1:nrow(pladb_fstat_alt_tab)

pladb_fdr_alt_tab <- pladb_fdr_alt[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(pladb_fdr_alt_tab) <- 1:nrow(pladb_fdr_alt_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(pladb_pdiff_alt_tab$EVENT, pladb_fstat_alt_tab$EVENT, pladb_fdr_alt_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_alt_diff_tab <- pladb_pdiff_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fstat_tab <- pladb_fstat_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fdr_tab <- pladb_fdr_alt_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_alt_diff_tab <- filtered_alt_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_alt_fstat_tab <- filtered_alt_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_alt_fdr_tab <- filtered_alt_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
pladb_alt_combination_tab <- filtered_alt_diff_tab %>%
  inner_join(filtered_alt_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_alt_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Render DataTable with enhancements
datatable(
  pladb_alt_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(pladb_alt_combination_tab,
          file = paste0("pladb_alt_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.3.3.2 PladB 1mM vs Control

pladb_pdiff_alt_tab <- pladb_1mm_control_pdiff_alt[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_pdiff_alt_tab) <- 1:nrow(pladb_pdiff_alt_tab)

pladb_fstat_alt_tab <- pladb_1mm_control_fstat_alt[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_fstat_alt_tab) <- 1:nrow(pladb_fstat_alt_tab)

pladb_fdr_alt_tab <- pladb_1mm_control_fdr_alt[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(pladb_fdr_alt_tab) <- 1:nrow(pladb_fdr_alt_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(pladb_pdiff_alt_tab$EVENT, pladb_fstat_alt_tab$EVENT, pladb_fdr_alt_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_alt_diff_tab <- pladb_pdiff_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fstat_tab <- pladb_fstat_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fdr_tab <- pladb_fdr_alt_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_alt_diff_tab <- filtered_alt_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_alt_fstat_tab <- filtered_alt_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_alt_fdr_tab <- filtered_alt_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
pladb_1mm_control_alt_combination_tab <- filtered_alt_diff_tab %>%
  inner_join(filtered_alt_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_alt_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Find if events appear in the control to 10mM pairwise comparison
pladb_1mm_control_alt_combination_tab$appears_in_10mm_control <- pladb_1mm_control_alt_combination_tab$EVENT %in% pladb_alt_combination_tab$EVENT

# Render DataTable with enhancements
datatable(
  pladb_1mm_control_alt_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(pladb_1mm_control_alt_combination_tab,
          file = paste0("pladb_1mm_control_alt_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.3.3.3 PladB 10mM vs 1mM

pladb_pdiff_alt_tab <- pladb_1mm_10mm_pdiff_alt[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_pdiff_alt_tab) <- 1:nrow(pladb_pdiff_alt_tab)

pladb_fstat_alt_tab <- pladb_1mm_10mm_fstat_alt[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(pladb_fstat_alt_tab) <- 1:nrow(pladb_fstat_alt_tab)

pladb_fdr_alt_tab <- pladb_1mm_10mm_fdr_alt[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(pladb_fdr_alt_tab) <- 1:nrow(pladb_fdr_alt_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(pladb_pdiff_alt_tab$EVENT, pladb_fstat_alt_tab$EVENT, pladb_fdr_alt_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_alt_diff_tab <- pladb_pdiff_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fstat_tab <- pladb_fstat_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fdr_tab <- pladb_fdr_alt_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_alt_diff_tab <- filtered_alt_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_alt_fstat_tab <- filtered_alt_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_alt_fdr_tab <- filtered_alt_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
pladb_10mm_1mm_alt_combination_tab <- filtered_alt_diff_tab %>%
  inner_join(filtered_alt_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_alt_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

pladb_10mm_1mm_alt_combination_tab$appears_in_10mm_control <- pladb_10mm_1mm_alt_combination_tab$EVENT %in% pladb_alt_combination_tab$EVENT

# Render DataTable with enhancements
datatable(
  pladb_10mm_1mm_alt_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(pladb_10mm_1mm_alt_combination_tab,
          file = paste0("pladb_10mm_1mm_alt_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.3.3.4 FMN2DKO

fmndko_pdiff_alt_tab <- fmndko_pdiff_alt[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(fmndko_pdiff_alt_tab) <- 1:nrow(fmndko_pdiff_alt_tab)

fmndko_fstat_alt_tab <- fmndko_fstat_alt[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(fmndko_fstat_alt_tab) <- 1:nrow(fmndko_fstat_alt_tab)

fmndko_fdr_alt_tab <- fmndko_fdr_alt[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(fmndko_fdr_alt_tab) <- 1:nrow(fmndko_fdr_alt_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(fmndko_pdiff_alt_tab$EVENT, fmndko_fstat_alt_tab$EVENT, fmndko_fdr_alt_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_alt_diff_tab <- fmndko_pdiff_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fstat_tab <- fmndko_fstat_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fdr_tab <- fmndko_fdr_alt_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_alt_diff_tab <- filtered_alt_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_alt_fstat_tab <- filtered_alt_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_alt_fdr_tab <- filtered_alt_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
fmndko_alt_combination_tab <- filtered_alt_diff_tab %>%
  inner_join(filtered_alt_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_alt_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Render DataTable with enhancements
datatable(
  fmndko_alt_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(fmndko_alt_combination_tab,
          file = paste0("fmndko_alt_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.3.3.5 Spire2

spire_pdiff_alt_tab <- spire_pdiff_alt[, c("GENE", "EVENT", "COORD", "Pdiff", "deltapsi")] %>%
  filter(!is.na(EVENT), Pdiff >= 0.90) %>%  # Remove rows where EVENT is NA and Pdiff >= 0.95
  arrange(desc(Pdiff), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(spire_pdiff_alt_tab) <- 1:nrow(spire_pdiff_alt_tab)

spire_fstat_alt_tab <- spire_fstat_alt[, c("GENE","EVENT","COORD","Fstat","deltapsi")] %>%
  filter(!is.na(EVENT)) %>%  # Remove rows where EVENT is NA
  arrange(desc(Fstat), desc(abs(deltapsi)))  # Sort by Pdiff and abs(deltapsi)
rownames(spire_fstat_alt_tab) <- 1:nrow(spire_fstat_alt_tab)

spire_fdr_alt_tab <- spire_fdr_alt[, c("GENE","EVENT","COORD","FDR","deltapsi")] %>%
  filter(!is.na(EVENT), FDR<=0.1) %>%  # Remove rows where EVENT is NA
  arrange((FDR), desc(abs(deltapsi)))
rownames(spire_fdr_alt_tab) <- 1:nrow(spire_fdr_alt_tab)

# Find common EVENTs across all tables
common_events <- Reduce(intersect, list(spire_pdiff_alt_tab$EVENT, spire_fstat_alt_tab$EVENT, spire_fdr_alt_tab$EVENT))

# Filter each table to include only rows with common EVENTs
filtered_alt_diff_tab <- spire_pdiff_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fstat_tab <- spire_fstat_alt_tab %>% filter(EVENT %in% common_events)
filtered_alt_fdr_tab <- spire_fdr_alt_tab %>% filter(EVENT %in% common_events)


# Select only unique columns from each table
filtered_alt_diff_tab <- filtered_alt_diff_tab %>% select(GENE,EVENT,COORD, deltapsi,Pdiff) # Extra column is Pdiff
filtered_alt_fstat_tab <- filtered_alt_fstat_tab %>% select(EVENT, Fstat) # Extra column is Fstat
filtered_alt_fdr_tab <- filtered_alt_fdr_tab %>% select(EVENT, FDR) # Extra column is FDR

# Merge the data frames by EVENT
spire_alt_combination_tab <- filtered_alt_diff_tab %>%
  inner_join(filtered_alt_fstat_tab, by = "EVENT") %>%
  inner_join(filtered_alt_fdr_tab, by = "EVENT") %>%
  distinct() %>%
  arrange(desc(abs(deltapsi)), desc(Fstat),desc(FDR), desc(Pdiff))

# Render DataTable with enhancements
datatable(
  spire_alt_combination_tab,
  options = list(
    pageLength = 10,                      # Rows per page
    autoWidth = TRUE,                     # Adjust column widths automatically
    dom = 'Bfrtip',                       # Add buttons for export
    buttons = c("copy", "csv", "excel", "pdf"),  # Simplified button definitions
    columnDefs = list(
      list(targets = "_all", className = "dt-center"),    # Center-align all columns
      list(
        targets = 3,  # Highlight significant Pdiff
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && parseFloat(data) > 0.99) {
               return '<span style=\"color: green; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      ),
      list(
        targets = 4,  # Highlight large deltapsi values
        render = JS(
          "function(data, type, row) {
             if (type === 'display' && Math.abs(parseFloat(data)) > 0.5) {
               return '<span style=\"color: red; font-weight: bold;\">' + data + '</span>';
             }
             return data;
           }"
        )
      )
    )
  ),
  rownames = FALSE,                       # Disable row names
  extensions = "Buttons"                  # Enable export buttons
) %>%
  formatStyle(
    columns = "Pdiff",
    color = styleInterval(0.99, c("black", "green")),
    fontWeight = styleEqual(0.99, "bold")  # Ensure consistent styling
  ) %>%
  formatStyle(
  columns = "deltapsi",
  color = styleInterval(c(-0.5, 0.5), c("black", "red", "black")),  # Three colors for two intervals
  fontWeight = styleInterval(c(-0.5, 0.5), c("normal", "bold", "normal"))  # Three weights for two intervals
)

write.csv(spire_alt_combination_tab,
          file = paste0("spire_alt_combination_table_", format(Sys.Date(), "%Y%m%d"), ".csv"),
          row.names = FALSE)

3.4 PSI Distribution Plots

These plots aim to explore any enrichment of either (intron/exon) retention or (intron/exon) exclusion. Symmetric plots show that inclusion and retention are equally common between conditions. Only showing events with |dPSI|>=0.1 for higher clarity. Transparency shows different studies.

3.4.1 PladB

# List of dataframes with metadata
dfs <- list(
  list(df = pladb_pdiff_exons, group = "Exons", study = "PladB_10mM_vs_Control"),
  list(df = pladb_1mm_control_pdiff_exons, group = "Exons", study = "PladB_1mM_vs_Control"),
  list(df = pladb_1mm_10mm_pdiff_exons, group = "Exons", study = "PladB_10mM_vs_1mm"),
  list(df = pladb_pdiff_introns, group = "Introns", study = "PladB_10mM_vs_Control"),
  list(df = pladb_1mm_control_pdiff_introns, group = "Introns", study = "PladB_1mM_vs_Control"),
  list(df = pladb_1mm_10mm_pdiff_introns, group = "Introns", study = "PladB_10mM_vs_1mm"),
   list(df = pladb_pdiff_alt, group = "Alt5/Alt3", study = "PladB_10mM_vs_Control"),
  list(df = pladb_1mm_control_pdiff_alt, group = "Alt5/Alt3", study = "PladB_1mM_vs_Control"),
  list(df = pladb_1mm_10mm_pdiff_alt, group = "Alt5/Alt3", study = "PladB_10mM_vs_1mm")
)

# Add metadata columns to each dataframe beforehand
dfs_with_metadata <- lapply(dfs, function(x) {
  x$df %>%
    mutate(Group = x$group, Study = x$study)
})

# Combine all dataframes into one using bind_rows() for faster merging
final_df <- bind_rows(dfs_with_metadata) %>%
  select(EVENT, Group, Study, deltapsi) %>%
  filter(abs(deltapsi)>=0.1) %>%
  na.omit()# Adjust this line if other columns should be prioritized


plot_psi_distribution <- ggplot(final_df, aes(x = deltapsi, fill = ifelse(deltapsi < 0, "Skipped", "Included"))) +
  geom_histogram( bins = 100, position = "identity") +
  labs(
    title = "PSI Distribution in PladB treated samples",
    x = "ΔPSI",
    y = "Density",
    fill = "|ΔPSI| >= 0.1",
    caption = paste0("Created by AG on ", Sys.Date())
  ) +
  theme_minimal(base_family = font) +
  theme(
    legend.position = "right",
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.text.x = element_text(size = 15),
    panel.background = element_rect(fill = "white"),
    panel.grid.major = element_line(color = "gray80", size = 0.5),
    panel.grid.minor = element_line(color = "gray95", size = 0.3),
    strip.text = element_text(size = 14, face = "bold"),
    panel.spacing = unit(1.5, "lines")  # Increase the distance between facets
  ) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~Group*Study)

plot_psi_distribution

3.4.2 FMN2DKO and SPIRE

# List of dataframes with metadata
dfs <- list(
  list(df = spire_pdiff_exons, group = "Exons", study = "Spire"),
  list(df = spire_pdiff_introns, group = "Introns", study = "Spire"),
  list(df = spire_pdiff_alt, group = "Alt5/Alt3", study = "Spire"),
  list(df = fmndko_pdiff_exons, group = "Exons", study = "FMN2DKO"),
  list(df = fmndko_pdiff_introns, group = "Introns", study = "FMN2DKO"),
  list(df = fmndko_pdiff_alt, group = "Alt5/Alt3", study = "FMN2DKO")
)

# Add metadata columns to each dataframe beforehand
dfs_with_metadata <- lapply(dfs, function(x) {
  x$df %>%
    mutate(Group = x$group, Study = x$study)
})

# Combine all dataframes into one using bind_rows() for faster merging
final_df <- bind_rows(dfs_with_metadata) %>%
  select(EVENT, Group, Study, deltapsi) %>%
  filter(abs(deltapsi)>=0.1) %>%
  na.omit()# Adjust this line if other columns should be prioritized


plot_psi_distribution <- ggplot(final_df, aes(x = deltapsi, fill = ifelse(deltapsi < 0, "Skipped", "Included"))) +
  geom_histogram( bins = 100, position = "identity") +
  labs(
    title = "PSI Distribution",
    x = "ΔPSI (DKO - WT)",
    y = "Density",
    fill = "|ΔPSI| >= 0.1",
    caption = paste0("Created by AG on ", Sys.Date())
  ) +
  theme_minimal(base_family = font) +
  theme(
    legend.position = "right",
    axis.title.x = element_text(margin = margin(t = 10)),
    axis.text.x = element_text(size = 15),
    panel.background = element_rect(fill = "white"),
    panel.grid.major = element_line(color = "gray80", size = 0.5),
    panel.grid.minor = element_line(color = "gray95", size = 0.3),
    strip.text = element_text(size = 14, face = "bold"),
    panel.spacing = unit(1.5, "lines")  # Increase the distance between facets
  ) +
  scale_fill_brewer(palette = "Set1") +
  facet_wrap(~Study*Group)

plot_psi_distribution

3.5 Enrichment Analysis with EnrichR

3.5.1 PladB 10mM vs Control

enrichdata<-enrichr(unique(c(pladb_exon_combination_tab$GENE, pladb_intron_combination_tab$GENE, pladb_alt_combination_tab$GENE)), databases = c("GO_Biological_Process_2023","GO_Cellular_Component_2023","GO_Molecular_Function_2023"))

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
##   Querying GO_Cellular_Component_2023... Done.
##   Querying GO_Molecular_Function_2023... Done.
## Parsing results... Done.

3.5.1.1 GO Biological Process

plotEnrich(enrichdata[[1]])

3.5.1.2 GO Cellular Component

plotEnrich(enrichdata[[2]])

3.5.1.3 GO Molecular Function

plotEnrich(enrichdata[[3]])

3.5.2 PladB 1mM vs Control

enrichdata<-enrichr(unique(c(pladb_1mm_control_exon_combination_tab$GENE, pladb_1mm_control_intron_combination_tab$GENE, pladb_1mm_control_alt_combination_tab$GENE)), databases = c("GO_Biological_Process_2023","GO_Cellular_Component_2023","GO_Molecular_Function_2023"))

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
##   Querying GO_Cellular_Component_2023... Done.
##   Querying GO_Molecular_Function_2023... Done.
## Parsing results... Done.

3.5.2.1 GO Biological Process

plotEnrich(enrichdata[[1]])

3.5.2.2 GO Cellular Component

plotEnrich(enrichdata[[2]])

3.5.2.3 GO Molecular Function

plotEnrich(enrichdata[[3]])

3.5.3 PladB 10mM vs 1mM

enrichdata<-enrichr(unique(c(pladb_10mm_1mm_exon_combination_tab$GENE, pladb_10mm_1mm_intron_combination_tab$GENE, pladb_10mm_1mm_alt_combination_tab$GENE)), databases = c("GO_Biological_Process_2023","GO_Cellular_Component_2023","GO_Molecular_Function_2023"))

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
##   Querying GO_Cellular_Component_2023... Done.
##   Querying GO_Molecular_Function_2023... Done.
## Parsing results... Done.

3.5.3.1 GO Biological Process

plotEnrich(enrichdata[[1]])

3.5.3.2 GO Cellular Component

plotEnrich(enrichdata[[2]])

3.5.3.3 GO Molecular Function

plotEnrich(enrichdata[[3]])

3.5.4 FMN2DKO

enrichdata<-enrichr(unique(c(fmndko_exon_combination_tab$GENE, fmndko_intron_combination_tab$GENE, fmndko_alt_combination_tab$GENE)), databases = c("GO_Biological_Process_2023","GO_Cellular_Component_2023","GO_Molecular_Function_2023"))

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
##   Querying GO_Cellular_Component_2023... Done.
##   Querying GO_Molecular_Function_2023... Done.
## Parsing results... Done.

3.5.4.1 GO Biological Process

plotEnrich(enrichdata[[1]])

3.5.4.2 GO Cellular Component

plotEnrich(enrichdata[[2]])

3.5.4.3 GO Molecular Function

plotEnrich(enrichdata[[3]])

3.5.5 Spire

enrichdata<-enrichr(unique(c(spire_exon_combination_tab$GENE, spire_intron_combination_tab$GENE, spire_alt_combination_tab$GENE)), databases = c("GO_Biological_Process_2023","GO_Cellular_Component_2023","GO_Molecular_Function_2023"))

## Uploading data to Enrichr... Done.
##   Querying GO_Biological_Process_2023... Done.
##   Querying GO_Cellular_Component_2023... Done.
##   Querying GO_Molecular_Function_2023... Done.
## Parsing results... Done.

3.5.5.1 GO Biological Process

plotEnrich(enrichdata[[1]])

3.5.5.2 GO Cellular Component

plotEnrich(enrichdata[[2]])

3.5.5.3 GO Molecular Function

plotEnrich(enrichdata[[3]])

4 Splicing Dynamics Heatmap PladB

set.seed(45)

# Extract all events
pladb_all <- filterEvents(
  pladb_events,
  types = c("C1", "C2", "C3", "S", "MIC", "IR","Alt3","Alt5"),
  N = 3
)

# Convert PSI values to a matrix and set row names
pladb_matrix <- pladb_all$PSI[, 7:15] %>%
  as.matrix()
rownames(pladb_matrix) <- pladb_all$PSI$EVENT

# Filter rows based on the final event lists and with abs(deltapsi)>=0.1
filtered_events <- unique(c(
  pladb_exon_combination_tab[abs(pladb_exon_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_1mm_control_exon_combination_tab[abs(pladb_1mm_control_exon_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_10mm_1mm_exon_combination_tab[abs(pladb_10mm_1mm_exon_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_intron_combination_tab[abs(pladb_intron_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_1mm_control_intron_combination_tab[abs(pladb_1mm_control_intron_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_10mm_1mm_intron_combination_tab[abs(pladb_10mm_1mm_intron_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_alt_combination_tab[abs(pladb_alt_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_1mm_control_alt_combination_tab[abs(pladb_1mm_control_alt_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_10mm_1mm_alt_combination_tab[abs(pladb_10mm_1mm_alt_combination_tab$deltapsi)>=0.1, "EVENT"]))
pladb_matrix <- pladb_matrix[rownames(pladb_matrix) %in% filtered_events, ]

# Arrange columns based on metadata
colnames(pladb_matrix)<-metadata_pladb$fastq_files
# Create a dendrogram for the columns
col_dend <- hclust(dist(t(pladb_matrix))) %>%
  as.dendrogram() %>%
  color_branches(k = 3) # Color branches with 3 clusters

# Check if row names contain "EX"
contains_EX <- grepl("EX", rownames(pladb_matrix))

# Define an improved color gradient for the heatmap
col_fun <- colorRamp2(
  c(0, 50, 100),
  c("purple", "white", "orange")
)

number_splices<-8


reference_table<-data.frame(EVENT=c(pladb_exon_combination_tab[abs(pladb_exon_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_1mm_control_exon_combination_tab[abs(pladb_1mm_control_exon_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_10mm_1mm_exon_combination_tab[abs(pladb_10mm_1mm_exon_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_intron_combination_tab[abs(pladb_intron_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_1mm_control_intron_combination_tab[abs(pladb_1mm_control_intron_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_10mm_1mm_intron_combination_tab[abs(pladb_10mm_1mm_intron_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_alt_combination_tab[abs(pladb_alt_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_1mm_control_alt_combination_tab[abs(pladb_1mm_control_alt_combination_tab$deltapsi)>=0.1, "EVENT"],
  pladb_10mm_1mm_alt_combination_tab[abs(pladb_10mm_1mm_alt_combination_tab$deltapsi)>=0.1, "EVENT"]), GENE=c(pladb_exon_combination_tab[abs(pladb_exon_combination_tab$deltapsi)>=0.1, "GENE"],
  pladb_1mm_control_exon_combination_tab[abs(pladb_1mm_control_exon_combination_tab$deltapsi)>=0.1, "GENE"],
  pladb_10mm_1mm_exon_combination_tab[abs(pladb_10mm_1mm_exon_combination_tab$deltapsi)>=0.1, "GENE"],
  pladb_intron_combination_tab[abs(pladb_intron_combination_tab$deltapsi)>=0.1, "GENE"],
  pladb_1mm_control_intron_combination_tab[abs(pladb_1mm_control_intron_combination_tab$deltapsi)>=0.1, "GENE"],
  pladb_10mm_1mm_intron_combination_tab[abs(pladb_10mm_1mm_intron_combination_tab$deltapsi)>=0.1, "GENE"],
  pladb_alt_combination_tab[abs(pladb_alt_combination_tab$deltapsi)>=0.1, "GENE"],
  pladb_1mm_control_alt_combination_tab[abs(pladb_1mm_control_alt_combination_tab$deltapsi)>=0.1, "GENE"],
  pladb_10mm_1mm_alt_combination_tab[abs(pladb_10mm_1mm_alt_combination_tab$deltapsi)>=0.1, "GENE"]))



ht<-Heatmap(
  pladb_matrix, name = "PSI",
  clustering_distance_rows = "pearson",
  col = col_fun,
  cluster_columns  = col_dend,
  column_title = "Oocytes PladB Samples at Different concentrations of PladB",
  row_title = "Splicing Events of |dPSI|>=0.1",
  column_names_rot = 45,
  column_labels = metadata_pladb$Description,
  column_dend_reorder = c(1:9),
  row_km = number_splices,
  row_dend_reorder = T,
  rect_gp = gpar(col = "white", lwd = 0.3),
  column_names_gp = gpar(fontsize = 16),
  show_row_names = FALSE,

)
ht = draw(ht); clustered_events<-row_order(ht)

text_list<-list()
# Loop through clustered events and assign dynamic names
for (i in 1:length(clustered_events)) {
  enrich_events <- enrichr(
    reference_table$GENE[reference_table$EVENT %in% rownames(pladb_matrix)[clustered_events[[i]]]],
    databases = c("GO_Biological_Process_2023", "GO_Cellular_Component_2023", "GO_Molecular_Function_2023")
  )

  # Assign the dynamic name and value
  text_list[[paste0("text", i)]] <- paste(
    enrich_events$GO_Cellular_Component_2023$Term[1],
    enrich_events$GO_Cellular_Component_2023$Term[2],
    enrich_events$GO_Cellular_Component_2023$Term[3],
    enrich_events$GO_Molecular_Function_2023$Term[1],
    enrich_events$GO_Molecular_Function_2023$Term[2],
    enrich_events$GO_Molecular_Function_2023$Term[3],
    sep = "; \n"
  )
}

protein_impact <- read.table("PROT_IMPACT-mm10-v2.3.tab", sep = "\t", 
                             header = TRUE, stringsAsFactors = FALSE, fill = TRUE, quote = "")

events_for_impact<-rownames(pladb_matrix[unlist(clustered_events, use.names = F),])

final_impact <- protein_impact %>%
  filter(EventID %in% events_for_impact) %>%
  arrange(match(EventID, events_for_impact)) %>%
  mutate(ONTO = ifelse(grepl("isoform", ONTO, ignore.case = TRUE), "1", ONTO)) %>%
  mutate(ONTO = ifelse(grepl("UTR", ONTO, ignore.case = TRUE), "-1", ONTO)) %>%
  mutate(ONTO = case_when(
    grepl("ORF", ONTO, ignore.case = TRUE) & grepl("inclusion", ONTO, ignore.case = TRUE) ~ "2",
    grepl("ORF", ONTO, ignore.case = TRUE) & grepl("exclusion", ONTO, ignore.case = TRUE) ~ "-2",
    TRUE ~ ONTO
  )) %>%
  mutate(ONTO = as.numeric(ONTO)) %>%
  replace_na(list(ONTO = 0))



# Define the left annotation (points) with axis labels
left_annotation <- rowAnnotation(
  Impact = anno_points(
    final_impact$ONTO,
    width = unit(4, "cm"),  # Adjusted width
    size = unit(4, "mm"),   # Adjusted size of points
    axis_param = list(
      side = "top",
      at = c(2, 1, 0, -1, -2), # Numeric values to label
      labels = c(
        "ORF disruption upon inclusion",
        "Alternative Isoform",
        "Unknown/NonCoding",
        "Regulatory",
        "ORF disruption upon exclusion"
      ),
      labels_rot = 45, # Rotation of labels
      gp = gpar(fontsize = 10) # Text style
    )
  )
)

# Define the right annotation (foo)
right_annotation <- rowAnnotation(
  foo = anno_empty(
    border = FALSE,
    width = max_text_width(unlist(text_list)) + unit(4, "mm")
  )
)

4.0.1 Heatmap showcasing the PSI of every candidate EX and IN

Note: Gene ontology of Cellular Component and Biological Process for each of the clustered blocks (splicing events that the pearson algorithm has found they behave similarly) appears on the right. The predicted impact on the protein appears on the left.

set.seed(45)

ht <- Heatmap(
  pladb_matrix, 
  name = "PSI",
  clustering_distance_rows = "pearson",
  col = col_fun,
  cluster_columns = col_dend,
  column_title = "Oocytes PladB Samples at Different Concentrations of PladB",
  column_title_gp = gpar(fontsize = 18, fontface = "bold"),
  row_title = "Splicing Events of |dPSI| >= 0.1",
  row_title_gp = gpar(fontsize = 16, fontface = "italic"),
  column_names_rot = 45,
  column_labels = metadata_pladb$Description,
  column_dend_reorder = c(1:9),
  row_km = number_splices,
  left_annotation = left_annotation,   # Add the points on the left
  right_annotation = right_annotation, # Add foo on the right
  row_dend_reorder = TRUE,
  column_names_gp = gpar(fontsize = 18, fontface = "bold", col = "darkblue"),
  show_row_names = FALSE,
  row_gap = unit(6, "mm"), # Adjust the gap size between rows
  column_gap = unit(4, "mm"), # Adjust the gap size between columns
  heatmap_legend_param = list(
    title = "PSI Value",
    title_gp = gpar(fontsize = 12),
    labels_gp = gpar(fontsize = 10),
    legend_position = c(-20, 0),  # Move the legend closer to the center (x, y coordinates)
    legend_direction = "vertical" # Align the legend horizontally
  )
)

# Draw the heatmap
ht <- draw(ht)

# Decorate the "foo" annotation slices (on the right)
for (i in seq_along(clustered_events)) {
  decorate_annotation("foo", slice = i, {
    grid.rect(x = 0, width = unit(2, "mm"), gp = gpar(fill = i, col = NA), just = "left")
    grid.text(paste(text_list[[i]], collapse = "\n"), x = unit(4, "mm"), just = "left", gp = gpar(fontsize = 17))
  })
}

5 System Settings

sessionInfo()

## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] grid      stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] dendextend_1.19.0      colorRamp2_0.1.0       ComplexHeatmap_2.21.1 
##  [4] clusterProfiler_4.10.1 DOSE_3.28.2            enrichplot_1.22.0     
##  [7] org.Hs.eg.db_3.18.0    org.Mm.eg.db_3.18.0    AnnotationDbi_1.64.1  
## [10] IRanges_2.36.0         S4Vectors_0.40.2       Biobase_2.62.0        
## [13] BiocGenerics_0.48.1    patchwork_1.3.0        enrichR_3.2           
## [16] biomaRt_2.58.2         readxl_1.4.3           ggtext_0.1.2          
## [19] showtext_0.9-7         showtextdb_3.0         sysfonts_0.8.9        
## [22] ggupset_0.4.0          paletteer_1.6.0        DT_0.33               
## [25] cowplot_1.1.3          lubridate_1.9.3        forcats_1.0.0         
## [28] stringr_1.5.1          purrr_1.0.2            readr_2.1.5           
## [31] tidyr_1.3.1            tibble_3.2.1           tidyverse_2.0.0       
## [34] dplyr_1.1.4            plotly_4.10.4          ggplot2_3.5.1         
## [37] betAS_1.2.1           
## 
## loaded via a namespace (and not attached):
##   [1] splines_4.3.1           later_1.4.1             prismatic_1.1.2        
##   [4] bitops_1.0-9            ggplotify_0.1.2         filelock_1.0.3         
##   [7] cellranger_1.1.0        polyclip_1.10-7         xts_0.14.1             
##  [10] XML_3.99-0.17           lifecycle_1.0.4         doParallel_1.0.17      
##  [13] lattice_0.22-6          MASS_7.3-60.0.1         crosstalk_1.2.1        
##  [16] backports_1.5.0         magrittr_2.0.3          sass_0.4.9             
##  [19] rmarkdown_2.29          jquerylib_0.1.4         yaml_2.3.10            
##  [22] rlist_0.4.6.2           httpuv_1.6.15           DBI_1.2.3              
##  [25] RColorBrewer_1.1-3      zlibbioc_1.48.2         ggraph_2.2.1           
##  [28] RCurl_1.98-1.16         yulab.utils_0.1.8       WriteXLS_6.7.0         
##  [31] tweenr_2.0.3            rappdirs_0.3.3          circlize_0.4.16        
##  [34] GenomeInfoDbData_1.2.11 ggrepel_0.9.6           tidytree_0.4.6         
##  [37] commonmark_1.9.2        highcharter_0.9.4       codetools_0.2-20       
##  [40] xml2_1.3.6              ggforce_0.4.2           shape_1.4.6.1          
##  [43] tidyselect_1.2.1        aplot_0.2.3             farver_2.1.2           
##  [46] viridis_0.6.5           matrixStats_1.4.1       BiocFileCache_2.10.2   
##  [49] jsonlite_1.8.9          GetoptLong_1.0.5        tidygraph_1.3.1        
##  [52] iterators_1.0.14        ggridges_0.5.6          foreach_1.5.2          
##  [55] tools_4.3.1             progress_1.2.3          treeio_1.31.0          
##  [58] Rcpp_1.0.13-1           glue_1.8.0              gridExtra_2.3          
##  [61] xfun_0.49               qvalue_2.34.0           TTR_0.24.4             
##  [64] GenomeInfoDb_1.38.8     withr_3.0.2             fastmap_1.2.0          
##  [67] fansi_1.0.6             digest_0.6.37           timechange_0.3.0       
##  [70] R6_2.5.1                mime_0.12               gridGraphics_0.5-1     
##  [73] colorspace_2.1-1        GO.db_3.18.0            markdown_1.13          
##  [76] RSQLite_2.3.9           utf8_1.2.4              generics_0.1.3         
##  [79] data.table_1.16.4       prettyunits_1.2.0       graphlayouts_1.2.1     
##  [82] httr_1.4.7              htmlwidgets_1.6.4       scatterpie_0.2.4       
##  [85] pkgconfig_2.0.3         gtable_0.3.6            blob_1.2.4             
##  [88] XVector_0.42.0          shadowtext_0.1.4        htmltools_0.5.8.1      
##  [91] fgsea_1.28.0            clue_0.3-66             scales_1.3.0           
##  [94] png_0.1-8               ggfun_0.1.8             knitr_1.49             
##  [97] rstudioapi_0.17.1       tzdb_0.4.0              reshape2_1.4.4         
## [100] rjson_0.2.23            nlme_3.1-166            curl_6.0.1             
## [103] GlobalOptions_0.1.2     cachem_1.1.0            zoo_1.8-12             
## [106] parallel_4.3.1          miniUI_0.1.1.1          HDO.db_0.99.1          
## [109] shinycssloaders_1.1.0   pillar_1.9.0            vctrs_0.6.5            
## [112] promises_1.3.2          dbplyr_2.5.0            cluster_2.1.4          
## [115] xtable_1.8-4            evaluate_1.0.1          cli_3.6.3              
## [118] compiler_4.3.1          rlang_1.1.4             crayon_1.5.3           
## [121] labeling_0.4.3          rematch2_2.1.2          plyr_1.8.9             
## [124] fs_1.6.5                stringi_1.8.4           viridisLite_0.4.2      
## [127] BiocParallel_1.36.0     assertthat_0.2.1        munsell_0.5.1          
## [130] Biostrings_2.70.3       lazyeval_0.2.2          colourpicker_1.3.0     
## [133] GOSemSim_2.28.1         Matrix_1.6-5            hms_1.1.3              
## [136] bit64_4.5.2             KEGGREST_1.42.0         shiny_1.9.1            
## [139] gridtext_0.1.5          igraph_2.1.1            broom_1.0.7            
## [142] memoise_2.0.1           bslib_0.8.0             thematic_0.1.6         
## [145] quantmod_0.4.26         ggtree_3.10.1           fastmatch_1.1-4        
## [148] bit_4.5.0.1             gson_0.1.0              ape_5.8